library(tidyverse)
library(ggpubr)
library(here)
library(patchwork)
temperature_range_data <- read_csv(here("1-data/TempChoice_abundance_dryad.csv"))
Temperature_ranges_plot <- ggplot(data = temperature_range_data, aes(x = day, y = mean_density, group = population, color = as.factor(temperature))) +
geom_point() +
geom_line() +
scale_y_log10() +
xlab("Day") +
ylab("Population density (cells/ml)") +
scale_color_manual(name = "Temperature (°C)",
values = c('#053061', '#2166ac', '#4393c3', '#92c5de',
'#d1e5f0','#fddbc7','#f4a582','#d6604d','#b2182b', '#67001f')) +
scale_x_continuous(breaks = c(5, 10)) +
theme(text = element_text(size = 15),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
legend.position = "bottom",
plot.margin = unit(c(1,1,1.5,1.2),"cm"))
Temperature_ranges_plot
# load the abundance data
abundance_data <- read_csv(here("1-data/TempAdapt_abundances_dryad.csv"))
# load the morphology data
morphology_data <- read_csv(here("1-data/TempAdapt_morphology_dryad.csv"))
# total number of cells measured
unique(morphology_data$id) %>% length()
## [1] 94344
# mean number of cells measured per day plus standard error
cells_summary <- morphology_data %>% group_by(population, date) %>% summarise(measured_cells = n())
cells_summary %>% ungroup() %>% summarise(mean_measured_cells = mean(measured_cells), se_measured_cells = sd(measured_cells)/sqrt(n()))
## # A tibble: 1 x 2
## mean_measured_cells se_measured_cells
## <dbl> <dbl>
## 1 513. 26.9
# calculate mean, variance and sd of area and aspect ratio (AR) for each population for each day
morphology_data_ts <- morphology_data %>% group_by(population, lineage, batch, date, day, temperature) %>%
summarise(mean_area_day = mean(mean_area), sd_area_day = sd(mean_area),
se_area_day = sd(mean_area)/sqrt(n()), var_area_day = var(mean_area),
mean_AR_day = mean(mean_ar), sd_AR_day = sd(mean_ar),
se_AR_day = sd(mean_ar)/sqrt(n()))
# calculate coefficient of variation of cell area
morphology_data_ts <- morphology_data_ts %>% mutate(cv_area_day = sd_area_day/mean_area_day)
# merge both datasets
morphology_abundance_data <- left_join(abundance_data, morphology_data_ts)
# add name column for plotting (lineage + temperature)
morphology_abundance_data <- morphology_abundance_data %>% mutate(name = paste0(lineage, "_", temperature))
# number of generations per population per batch
generations <- morphology_abundance_data %>% group_by(population, lineage, batch,replicate, temperature) %>%
summarise(maxi_density = max(mean_density), min_density = min(mean_density),
number_generations = (log(max(mean_density)) - log(min(mean_density)))/log(2))
# mean number of generations per batch per temperature
generations_summary <- generations %>% group_by(batch, temperature) %>% summarise(mean_generations_batch = mean(number_generations))
# create a tibble with the number of generations per batch at each temperature
generations_summary <- generations_summary %>% group_by(batch) %>% pivot_wider(names_from = temperature, values_from = mean_generations_batch, names_prefix = "generations_")
# round the number of generations
generations_summary <- generations_summary %>% mutate(text_plot = case_when(
batch == 1 ~ paste0(round(generations_20, digits = 0), " gen."),
batch == 2 | batch == 3 ~ paste0(round(generations_38, digits = 0), " gen."),
batch > 3 ~ paste0(round(generations_20, digits = 0), " gen. 20°C \n", round(generations_38, digits = 0), " gen. 38°C")))
# define size of small and large fonts for the plot
small_font <- 10
large_font <- 13
# plot population abundances through time
plot_abundances <- ggplot(data = morphology_abundance_data) +
geom_point(aes(x = batch_day, y = mean_density, group = population, colour = name), size = 1) +
geom_line(aes(x = batch_day, y = mean_density, group = population, colour = name)) +
geom_label(data = generations_summary, aes(x = 9, y = 500, label = text_plot), size = 2.1) +
scale_y_log10() +
ylab("Population density\n(cells/ml)") +
xlab("") +
scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b", "red", "blue"),
guide = F, name = "Population and\n temperature",
labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
"2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
"3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
"4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
scale_x_continuous(breaks = c(5, 10)) +
theme(text = element_text(size = large_font),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
legend.position = "bottom",
plot.margin = unit(c(0, 1.5, 0 , 1.5 ),"cm")) +
facet_wrap(~ batch, ncol = 5, labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2",
"3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5")))
# minimum and maximum cell area in the control
area_min <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>%
select(mean_area_day) %>% min()
area_max <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>%
select(mean_area_day) %>% max()
# plot cell area
plot_cell_area <- ggplot(data = morphology_abundance_data, aes(x = batch_day, y = mean_area_day,
group = population, col = name)) +
geom_point(size = 1) +
geom_hline(aes(yintercept = area_min), linetype = "dashed") +
geom_hline(aes(yintercept = area_max), linetype = "dashed") +
geom_errorbar(aes(ymin = mean_area_day - se_area_day, ymax = mean_area_day + se_area_day, width = 0.9), size = 0.2) +
geom_line() +
xlab("") +
ylab("Cell size (um2)") +
scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
name = "Population and\n temperature",
labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
"2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
"3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
"4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
scale_x_continuous(breaks = c(5, 10)) +
theme(text = element_text(size = large_font), legend.text = element_text(size = small_font),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
legend.position = "bottom", strip.text = element_blank(),
plot.margin = unit(c(0, 1.5, 0 ,1.5 ),"cm")) +
facet_wrap(~ batch, ncol = 5, labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2",
"3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5")))
# minimum and maximum cv of cell area in the control
cv_area_min <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>%
select(cv_area_day) %>% min()
cv_area_max <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>%
select(cv_area_day) %>% max()
# plot coefficient of variation of cell area
# multiply cv per 100 to obtain percentage
plot_cv_cell_area <- ggplot(data = morphology_abundance_data, aes(x = batch_day, y = cv_area_day*100,
group = population, col = name)) +
geom_point(size = 1) +
geom_hline(aes(yintercept = cv_area_min*100), linetype = "dashed") +
geom_hline(aes(yintercept = cv_area_max*100), linetype = "dashed") +
geom_line() +
xlab("") +
ylab("Coefficient of variation\n of cell size (%)") +
scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
name = "Population and\n temperature",
labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
"2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
"3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
"4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
scale_x_continuous(breaks = c(5, 10)) +
scale_y_continuous(breaks = c(15, 30, 45, 60)) +
theme(text = element_text(size = large_font), legend.text = element_text(size = small_font),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
legend.position = "bottom", strip.text = element_blank(),
plot.margin = unit(c(0, 1.5, 0, 1.5 ),"cm")) +
facet_wrap(~ batch, ncol = 5, labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2",
"3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5")))
# minimum and maximum variance of cell area in the control
var_area_min <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>%
select(var_area_day) %>% min()
var_area_max <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>%
select(var_area_day) %>% max()
# plot variance of cell area
plot_var_cell_area <- ggplot(data = morphology_abundance_data, aes(x = batch_day, y = var_area_day, group = population, col = name)) +
geom_point(size = 1) +
geom_hline(aes(yintercept = var_area_min), linetype = "dashed") +
geom_hline(aes(yintercept = var_area_max), linetype = "dashed") +
geom_line() +
xlab("Day") +
ylab("Variance of\n cell size") +
scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
name = "Population and\n temperature",
labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
"2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
"3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
"4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
scale_x_continuous(breaks = c(5, 10)) +
theme(text = element_text(size = large_font), legend.text = element_text(size = small_font),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
legend.position = "bottom", strip.text = element_blank(),
plot.margin = unit(c(0, 1.5, 0 ,1.5 ),"cm")) +
facet_wrap(~ batch, ncol = 5, labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2",
"3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5")))
# minimum and maximum cell aspect ratio in the control
AR_min <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>%
select(mean_AR_day) %>% min()
AR_max <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>%
select(mean_AR_day) %>% max()
# plot cell aspect ratio
plot_cell_AR <- ggplot(data = morphology_abundance_data, aes(x = batch_day, y = mean_AR_day,
group = population, col = name)) +
geom_point(size = 1) +
geom_hline(aes(yintercept = AR_min), linetype = "dashed") +
geom_hline(aes(yintercept = AR_max), linetype = "dashed") +
geom_errorbar(aes(ymin = mean_AR_day - se_AR_day, ymax = mean_AR_day + se_AR_day, width = 0.9), size = 0.2) +
geom_line() +
xlab("") +
ylab("Cell shape") +
scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
name = "Population and\n temperature",
labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
"2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
"3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
"4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
scale_x_continuous(breaks = c(5, 10)) +
theme(text = element_text(size = large_font), legend.text = element_text(size = small_font),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
legend.position = "bottom", strip.text = element_blank(),
plot.margin = unit(c(0, 1.5, 0 ,1.5 ),"cm")) +
facet_wrap(~ batch, ncol = 5, labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2",
"3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5")))
# all plots in one figure
all_time_series <- ggarrange(plot_abundances, plot_cell_area, plot_var_cell_area, plot_cv_cell_area, plot_cell_AR,
nrow = 5, common.legend = T, legend = "bottom", align = "hv",
labels = "AUTO")
# print the plot
all_time_series
# load the package growthrates
library(growthrates)
# select only relevant data for growth model (population, population density and temperature)
# exclude populations 32, 141 and 341, since model will be run separetely for them
data_subset1 <- morphology_abundance_data %>%
group_by(population, lineage, batch, replicate, temperature) %>%
select(population, mean_density, batch_day) %>%
filter(!population %in% c("32", "141", "341"))
# set parameters to start the model (p), lower boundary and upper boundary
# y0 = initial population density, mumax = maximum growth rate, K = carrying capacity, lambda = lag phase
p <- c(y0 = 100, mumax = 1, K = 1000000, lambda = 5)
lower <- c(y0 = 0, mumax = 0, K = 100000, lambda = 0)
upper <- c(y0 = 5000, mumax = 15, K = 5000000, lambda = 8)
# run gompertz model with lag phase estimation for each population separately
subset1_gompertz3 <- all_growthmodels(
mean_density ~ grow_gompertz3(batch_day, parms) | population,
data = data_subset1,
p = p, lower = lower, upper = upper,
log = "y", ncores = 4)
# check the rsquared of the models
# r squared of the model
rsquared(subset1_gompertz3) > 0.9
rsquared(subset1_gompertz3) # all above 0.90, expect pop. 241, with 0.80
# coefficients of the model
coef(subset1_gompertz3)
# plot the model, y in log scale
par(mfrow = c(8, 4))
par(mar = c(2.5, 4, 2, 1))
plot(subset1_gompertz3)
plot(subset1_gompertz3, log = "y")
# populations 32, 141, 341 presented suboptimal fittings with the above parameters
# model these populations separately
# population 32, remove last day since population is already crashing
pop32 <- morphology_abundance_data %>%
group_by(population, lineage, batch, replicate, temperature) %>%
select(population, mean_density, batch_day) %>%
filter(population == 32 & batch_day < 12)
# same parameters as previous populations
pop32_gompertz <- fit_growthmodel(FUN = grow_gompertz3, p = p,
pop32$batch_day, pop32$mean_density,
lower = lower, upper = upper)
# # population 141
pop141 <- morphology_abundance_data %>%
group_by(population, lineage, batch, replicate, temperature) %>%
select(population, mean_density, batch_day) %>%
filter(population == 141 & batch_day < 6)
# use new parameters
p2 <- c(y0 = 500, mumax = 1, K = 1000000, lambda = 0.1)
lower2 <- c(y0 = 0, mumax = 0, K = 100000, lambda = 0)
upper2 <- c(y0 = 10000, mumax = 30, K = 5000000, lambda = 4)
# run model for population 141
pop141_gompertz <- fit_growthmodel(FUN = grow_gompertz3, p = p,
pop141$batch_day, pop141$mean_density,
lower = lower, upper = upper)
# population 341, remove final day that population is crashing
pop341 <- morphology_abundance_data %>%
group_by(population, lineage, batch, replicate, temperature) %>%
select(population, mean_density, batch_day) %>%
filter(population == 341 & batch_day <= 5)
# rsquared of the model is very low, manually estimate growth rate and lag phase
# select relevant days
pop341_day0 <- pop341 %>% filter(batch_day == 0)
pop341_day3 <- pop341 %>% filter(batch_day == 3)
# calculate maximum growth rate as the log difference in population density divided by the time difference
mumax <- (log10(pop341_day3$mean_density) - log10(pop341_day0$mean_density))/(pop341_day3$batch_day - pop341_day0$batch_day)
lambda <- 0 # define lag phase as 0
y0 <- NA
K <- NA
population <- 341
# store growth parameters in single object
pop341_manual_growth <- tibble(mumax, lambda, y0, K, population)
remove(mumax, lambda, y0, K, population)
# check model fit and coeficients
# pop 32
rsquared(pop32_gompertz)
coef(pop32_gompertz)
par(mfrow = c(1,2))
plot(pop32_gompertz)
plot(pop32_gompertz, log = "y")
# pop 141
rsquared(pop141_gompertz)
coef(pop141_gompertz)
par(mfrow = c(1,2))
plot(pop141_gompertz)
plot(pop141_gompertz, log = "y")
# get the model predictions for all gompertz models
# new time variable with the maximum batch_day across all populations (12)
time <- data.frame(time = seq(from = 0, to = 12, by = 1))
# predict function for each problematic population, add name of the population to the tibble
pop32_growth_prediction <- pop32_gompertz %>% predict(newdata = time) %>% as_tibble() %>% dplyr::rename(batch_day = time, mean_density_predict = y) %>% mutate(population = "32")
pop141_growth_prediction <- pop141_gompertz %>% predict(newdata = time) %>% as_tibble() %>% dplyr::rename(batch_day = time, mean_density_predict = y) %>% mutate(population = "141")
# predict function for each the remaining populations
subset1_growth_prediction <- predict(subset1_gompertz3, newdata = time)
# transform list into data frame with population as a column
subset1_growth_prediction <- map_df(subset1_growth_prediction, ~as.data.frame(.x), .id = "population")
# rename the columns
subset1_growth_prediction <- subset1_growth_prediction %>% dplyr::rename(batch_day = time, mean_density_predict = y)
# merge all models
prediction_gompertz <- bind_rows(subset1_growth_prediction, pop32_growth_prediction, pop141_growth_prediction)
prediction_gompertz$population <- as.numeric(prediction_gompertz$population)
# add growth model predictions to main dataset
morphology_abundance_data_gompertz <- morphology_abundance_data %>% left_join(prediction_gompertz)
# create a dataset with the coeficients of the model
coef_gompertz <- coef(subset1_gompertz3) %>% as_tibble()
coef_gompertz <- cbind(coef_gompertz, population = as.double(names(subset1_gompertz3)))
coef_gompertz <- coef_gompertz %>% rbind(c(coef(pop32_gompertz), 32))
coef_gompertz <- coef_gompertz %>% rbind(c(coef(pop141_gompertz), 141))
coef_gompertz <- coef_gompertz %>% rbind(pop341_manual_growth)
# select metadata from the original dataset
metadata <- morphology_abundance_data %>% select(population, lineage, batch, replicate, temperature) %>% distinct()
# add metadata to the growth model coeficients data
coef_gompertz <- coef_gompertz %>% left_join(metadata)
# add name column for plotting (lineage + temperature)
coef_gompertz <- coef_gompertz %>% mutate(name = paste0(lineage, "_", temperature))
# define font sizes
small_font <- 10
large_font <- 13
# plot the lag phase length (lambda)
plot_lambda <- ggplot(coef_gompertz, aes(x = batch, y = lambda, fill = name, group = population)) +
geom_col(position = "dodge") +
xlab("") +
ylab("Lag phase (days)") +
scale_fill_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
name = "Population and\n temperature",
labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
"2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
"3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
"4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
guides(fill = guide_legend(title.theme = element_text(size = small_font), label.theme = element_text(size = small_font))) +
theme(text = element_text(size = large_font),
panel.border = element_rect(colour = "black", fill = NA),
panel.background = element_blank(),
legend.position = "right",
#plot.margin = margin(t = 10, r = 10, b = 10, l = 10),
axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.text.x = element_blank(), axis.ticks.x = element_blank()) + #,
facet_grid(~batch, scales = "free_x",
labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2",
"3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5")))
# plot the maximum growth rate (mumax)
plot_mumax <- ggplot(coef_gompertz, aes(x = batch, y = mumax, fill = name, group = population)) +
geom_col(position = "dodge") +
xlab("") +
ylab(bquote("Max. growth rate"~(day^-1))) +
scale_fill_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
name = "Population and\n temperature", guide = F,
labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
"2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
"3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
"4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
#guides(fill = guide_legend(title.theme = element_text(size = small_font), label.theme = element_text(size = small_font))) +
theme(text = element_text(size = large_font),
panel.border = element_rect(colour = "black", fill = NA),
panel.background = element_blank(),
legend.position = "bottom",
#plot.margin = margin(t = 10, r = 10, b = 10, l = 10),
axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.text.x = element_blank(), axis.ticks.x = element_blank()) + #,
facet_grid(~batch, scales = "free_x",
labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2",
"3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5")))
gompertz_fit_plot <- ggplot(data = morphology_abundance_data_gompertz) +
geom_point(aes(x = batch_day, y = mean_density, group = population, colour = name), size = 1.5) +
geom_line(aes(x = batch_day, y = mean_density_predict, group = population, colour = name)) +
scale_y_log10() +
ylab("Population density (cells/ml)\n") +
xlab("Day\n") +
scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b", "red", "blue"),
name = "Population and\n temperature", guide = F,
labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
"2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
"3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
"4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
scale_x_continuous(breaks = c(5, 10)) +
theme(text = element_text(size = large_font),
panel.background = element_blank(), panel.border = element_rect(colour = "black", fill = NA),
axis.line = element_line(colour = "black"),
legend.position = "bottom",
plot.margin = unit(c(0, 1.5, 0 , 1.5 ),"cm")) +
facet_grid(paste0("Batch ", batch) ~ paste0("Population ", lineage))
# layout of the final figure
layout <- "
AA#
AA#
AA#
BBD
CCD"
# create final figure with three facets
plot_growth <- gompertz_fit_plot + plot_mumax + plot_lambda + guide_area() + plot_layout(guides = "collect", design = layout) + plot_annotation(tag_levels = 'A')
# print figure
plot_growth
# create variable merging batch, temperature and replicate; select only used variables
growth_response_gompertz <- coef_gompertz %>% ungroup() %>%
mutate(batch_temp = paste0("b", batch, "_", temperature, "_", replicate)) %>%
select(-population, -batch, - replicate, - temperature, -name)
# calculate % change in comparison to the same population lineage in batch 1
growth_response_gompertz <- growth_response_gompertz %>% select(-y0, -K) %>% pivot_longer(cols = 1:2, names_to = "data_type")
growth_response_gompertz <- growth_response_gompertz %>% group_by(lineage, data_type) %>% pivot_wider(names_from = batch_temp, values_from = value)
# calculate the difference in maximum growth rate and in lag phase duration in absolute values, not the % change since numbers are too large and interpretation is difficult
growth_response_gompertz <- growth_response_gompertz %>%
mutate(week2 = b2_38_NA - b1_20_NA,
week3_1 = b3_38_1 - b1_20_NA,
week3_2 = b3_38_2 - b1_20_NA,
week4 = b4_38_1 - b1_20_NA,
week5 = b5_38_1 - b1_20_NA,
week4_cold = b4_20_2 - b1_20_NA,
week5_cold = b5_20_2 - b1_20_NA)
# tidy the data
growth_response_gompertz <- growth_response_gompertz %>% select(1, 2, 11:17) %>%
gather(key = "comparison", value = "percent_batch1", 3:9)
# repeat the same calculations with the morphological traits
# calculate mean cell area and cell aspect ratio in each batch
morphology_data_batch <- morphology_abundance_data %>% drop_na(mean_area_day) %>%
group_by(population, lineage, batch, replicate, temperature) %>%
summarise(mean_area_batch = mean(mean_area_day), cv_area_batch = mean(cv_area_day), var_area_batch = mean(var_area_day),
mean_AR_batch = mean(mean_AR_day))
# create variable combining batch, temperature and replicate; select only used variables
morphology_response <- morphology_data_batch %>% ungroup() %>%
mutate(batch_temp = paste0("b",batch, "_", temperature, "_", replicate)) %>%
select(lineage, batch_temp, mean_area_batch, mean_AR_batch, cv_area_batch, var_area_batch)
morphology_response <- morphology_response %>% gather(key = "data_type", value = "value", 3:6)
morphology_response <- morphology_response %>% group_by(lineage, data_type) %>% spread(key = batch_temp, value = value)
# calculate % change in comparison to the same population lineage in batch 1
morphology_response <- morphology_response %>% group_by(lineage, data_type) %>%
mutate(week2 = ((b2_38_NA - b1_20_NA) * 100) / b1_20_NA,
week3_1 = ((b3_38_1 - b1_20_NA) * 100) / b1_20_NA,
week3_2 = ((b3_38_2 - b1_20_NA) * 100) / b1_20_NA,
week4 = ((b4_38_1 - b1_20_NA) * 100) / b1_20_NA,
week5 = ((b5_38_1 - b1_20_NA) * 100) / b1_20_NA,
week4_cold = ((b4_20_2 - b1_20_NA) * 100) / b1_20_NA,
week5_cold = ((b5_20_2 - b1_20_NA) * 100) / b1_20_NA)
morphology_response <- morphology_response %>% select(1, 2, 11:17) %>%
gather(key = "comparison", value = "percent_batch1", 3:9)
# merge growth and lag phase with the morphology data
response_data <- rbind.data.frame(growth_response_gompertz, morphology_response)
# add batch information
response_data <- response_data %>% ungroup() %>% mutate(batch = case_when(
comparison == "week2" ~ 2,
comparison == "week3_1" ~ 3,
comparison == "week3_2" ~ 3,
comparison == "week4" ~ 4,
comparison == "week5" ~ 5,
comparison == "week4_cold" ~ 4,
comparison == "week5_cold" ~ 5
))
# add temperature information
response_data <- response_data %>% mutate(temperature = case_when(
comparison == "week2" ~ 38,
comparison == "week3_1" ~ 38,
comparison == "week3_2" ~ 38,
comparison == "week4" ~ 38,
comparison == "week5" ~ 38,
comparison == "week4_cold" ~ 20,
comparison == "week5_cold" ~ 20
))
# add variable for lineage and temperature color
response_data <- response_data %>% mutate(name = paste0(lineage, "_", temperature))
# change data type into factor and change levels for a nicer looking plot
response_data$data_type <- factor(response_data$data_type, levels = c("mumax", "lambda",
"mean_area_batch", "cv_area_batch", "var_area_batch",
"mean_AR_batch"))
# transform lineage into categorical variable
response_data$lineage <- as.character(response_data$lineage)
# look at the data
ggplot(data = response_data, aes(x = batch, y = percent_batch1, color = name)) +
geom_point(alpha = 0.8, size = 2) +
geom_hline(yintercept = 0, linetype = "dashed") +
scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b", "#878787"),
name = "Population and\n temperature",
labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
"2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
"3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
"4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
ylab("Change compared to control (%)") +
xlab("Batch culture") +
guides(fill = guide_legend(title.theme = element_text(size = 15),
label.theme = element_text(size = 15))) +
theme(text = element_text(size = 15),
panel.border = element_rect(colour = "black", fill = NA),
panel.background = element_blank(),
legend.position = "bottom") +
facet_wrap(~data_type, ncol = 3, labeller = as_labeller(c("mumax" = "Growth rate",
"lambda" = "Lag phase",
"cv_area_batch" = "CV cell size",
"mean_AR_batch" = "Cell shape",
"mean_area_batch" = "Cell size",
"var_area_batch" = "Variance cell size")), scales = "free_y")
# load the libraries
library(MCMCglmm)
library(coda)
library(purrr)
# set seed
set.seed(1)
# select the populations that experienced 38 °C
response_data_38 <- response_data %>% filter(temperature == 38)
# create variable for time difference since warming started
response_data_38 <- response_data_38 %>% mutate(batches_38 = batch - 2)
# create a list with the prior information, using an uninformative prior
prior_1 <- list(R = list(V = 1, nu = 0.002),
G = list(G1 = list(V = 1, nu = 0.002)))
# run one separate model for each of the 5 response variables
models_38 <- response_data_38 %>% group_by(data_type) %>% nest() %>%
mutate(mixed_model = map(data, ~MCMCglmm(fixed = percent_batch1 ~ batches_38,
random = ~lineage,
data = .,
family = "gaussian",
prior = prior_1,
nitt = 2000000,
burnin = 30000,
thin = 1000, verbose = F)))
# create new data
new_data <- data.frame(expand.grid(lineage = c("1"), batches_38 = c(0:3),
percent_batch1 = 0))
# get confidence interval from each model
models_38 <- models_38 %>% mutate(predictions = map(mixed_model, ~predict.MCMCglmm(object = .,
newdata = new_data,
interval = "confidence")))
models_38 <- models_38 %>% mutate(predictions = map(.x = predictions, ~as_tibble(.)))
# merge new data and model fit for plotting
models_38 <- models_38 %>% mutate(predictions = map(.x = predictions, ~cbind(new_data, .)))
# store models in a separate object
model_38_predictions <- unnest(models_38, predictions)
# plot the data and the model
response_38_plot <- ggplot(data = response_data_38, aes(x = batch, y = percent_batch1, color = name)) +
geom_point(aes(color = as.factor(lineage)), alpha = 0.8, size = 2) +
geom_smooth(data = model_38_predictions, aes(x = batches_38 + 2, y = fit, ymin = lwr, ymax = upr), colour = "black", stat = "identity") +
geom_hline(yintercept = 0, linetype = "dashed") +
#scale_y_continuous(limits = c(-95, 95)) +
xlab("Batch") +
ylab("Change in comparison to control (%)") +
scale_color_manual(values = c("#efd0d4", "#d88b95","#c14655","#b2182b"),
name = "Population", labels = as_labeller(c( "1" = "Pop. 1",
"2" = "Pop. 2",
"3" = "Pop. 3",
"4" = "Pop. 4"))) +
theme(text = element_text(size = 15), plot.title = element_text(hjust = 0.5, margin = unit(c(0.5,0.5,0.5,0.5), "lines")),
panel.border = element_rect(colour = "black", fill = NA),
panel.background = element_blank(),
legend.position = "bottom") +
facet_wrap(~data_type, ncol = 3, labeller = as_labeller(c("mumax" = "Growth rate",
"lambda" = "Lag phase",
"cv_area_batch" = "CV cell size",
"mean_AR_batch" = "Cell shape",
"mean_area_batch" = "Cell size",
"var_area_batch" = "Variance cell size")), scales = "free_y")
# print the plot
response_38_plot
# select the populations that experienced 20 °C
response_data_20 <- response_data %>% filter(temperature == 20)
# create variable for time difference since return to 20 °C started
response_data_20 <- response_data_20 %>% mutate(batches_20 = batch - 4)
# create a list with the prior information, using an uninformative prior
# same prior as 38°C models
prior_1_new <- list(R = list(V = 1, nu = 0.002),
G = list(G1 = list(V = 1, nu = 0.002)))
# run one separate model for each of the 5 response variables
models_20 <- response_data_20 %>% group_by(data_type) %>% nest() %>%
mutate(mixed_model = map(data, ~MCMCglmm(fixed = percent_batch1 ~ batches_20,
random = ~lineage,
data = .,
family = "gaussian",
prior = prior_1_new,
nitt = 2000000,
burnin = 30000,
thin = 1000,
verbose = F)))
# create new data
new_data_20 <- data.frame(expand.grid(lineage = c("1"), batches_20 = c(0, 1),
percent_batch1 = 0))
# get confidence interval from each model
models_20 <- models_20 %>% mutate(predictions = map(mixed_model, ~predict.MCMCglmm(object = .,
newdata = new_data_20,
interval = "confidence")))
models_20 <- models_20 %>% mutate(predictions = map(.x = predictions, ~as_tibble(.)))
# merge new data and model fit for plotting
models_20 <- models_20 %>% mutate(predictions = map(.x = predictions, ~cbind(new_data_20, .)))
# store models in a separate object
model_20_predictions <- unnest(models_20, predictions)
# plot the data and the model
response_20_plot <- ggplot(data = response_data_20, aes(x = batch, y = percent_batch1, color = name)) +
geom_point(aes(color = as.factor(lineage)), alpha = 0.8, size = 2) +
geom_smooth(data = model_20_predictions,
aes(x = batches_20 + 4, y = fit, ymin = lwr, ymax = upr),
colour = "black", stat = "identity") +
geom_hline(yintercept = 0, linetype = "dashed") +
scale_x_continuous(breaks = c(4, 5)) +
xlab("Batch") +
ylab("Change in comparison to control (%)") +
scale_color_manual(values = c("#d2e0ee", "#90b2d5", "#4d84bc", "#2166ac"),
name = "Population",
labels = as_labeller(c( "1" = "Pop. 1",
"2" = "Pop. 2",
"3" = "Pop. 3",
"4" = "Pop. 4"))) +
theme(text = element_text(size = 15), plot.title = element_text(hjust = 0.5, margin = unit(c(0.5,0.5,0.5,0.5), "lines")),
panel.border = element_rect(colour = "black", fill = NA),
panel.background = element_blank(),
legend.position = "bottom") +
facet_wrap(~data_type, ncol = 3, labeller = as_labeller(c("mumax" = "Growth rate",
"lambda" = "Lag phase",
"cv_area_batch" = "CV cell size",
"mean_AR_batch" = "Cell shape",
"mean_area_batch" = "Cell size",
"var_area_batch" = "Variance cell size")), scales = "free_y")
# print the plot
response_20_plot
# store MCMC results in a separate object
models_38_MCMC <- models_38[[3]]
# name each MCMC run with the variable name
models_38_MCMC <- set_names(models_38_MCMC, models_38$data_type)
# check autocorrelation of the iterations, should be less than 0.1
map(models_38_MCMC, "Sol") %>% as.mcmc() %>% autocorr()
map(models_38_MCMC, "VCV") %>% as.mcmc() %>% autocorr()
# plot trace
plot(models_38_MCMC$mumax$Sol)
plot(models_38_MCMC$lambda$Sol)
plot(models_38_MCMC$mean_area_batch$Sol)
plot(models_38_MCMC$cv_area_batch$Sol)
plot(models_38_MCMC$var_area_batch$Sol)
plot(models_38_MCMC$mean_AR_batch$Sol)
plot(log10(models_38_MCMC$mumax$VCV))
plot(log10(models_38_MCMC$lambda$VCV))
plot(log10(models_38_MCMC$mean_area_batch$VCV))
plot(log10(models_38_MCMC$cv_area_batch$VCV))
plot(log10(models_38_MCMC$var_area_batch$VCV))
plot(log10(models_38_MCMC$mean_AR_batch$VCV))
# store MCMC results in a separate object
models_20_MCMC <- models_20[[3]]
# name each MCMC run with the variable name
models_20_MCMC <- set_names(models_20_MCMC, models_20$data_type)
# check autocorrelation of the iterations, should be less than 0.1
map(models_20_MCMC, "Sol") %>% as.mcmc() %>% autocorr()
map(models_20_MCMC, "VCV") %>% as.mcmc() %>% autocorr()
# plot trace
plot(models_20_MCMC$mumax$Sol)
plot(models_20_MCMC$lambda$Sol)
plot(models_20_MCMC$mean_area_batch$Sol)
plot(models_20_MCMC$cv_area_batch$Sol)
plot(models_20_MCMC$var_area_batch$Sol)
plot(models_20_MCMC$mean_AR_batch$Sol)
plot(log10(models_20_MCMC$mumax$VCV))
plot(log10(models_20_MCMC$lambda$VCV))
plot(log10(models_20_MCMC$mean_area_batch$VCV))
plot(log10(models_20_MCMC$cv_area_batch$VCV))
plot(log10(models_20_MCMC$var_area_batch$VCV))
plot(log10(models_20_MCMC$mean_AR_batch$VCV))
# plot both models in the same plot with all the response data
response_plot <- ggplot(data = filter(response_data, data_type != "K"), aes(x = batch, y = percent_batch1, color = name)) +
geom_point(alpha = 0.8, size = 2) +
geom_hline(yintercept = 0, linetype = "dashed") +
scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
name = "Population and\n temperature",
labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
"2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
"3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
"4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
geom_smooth(data = model_20_predictions %>% filter(data_type != "K"),
aes(x = batches_20 + 4, y = fit, ymin = lwr, ymax = upr), alpha = 0.2, colour = "#82a8d0", stat = "identity") +
geom_smooth(data = model_38_predictions %>% filter(data_type != "K"),
aes(x = batches_38 + 2, y = fit, ymin = lwr, ymax = upr), alpha = 0.2, colour = "#d37d88", stat = "identity") +
geom_hline(yintercept = 0, linetype = "dashed") +
xlab("Batch") +
ylab("Change in comparison to control (%)") +
theme(text = element_text(size = 15),
panel.border = element_rect(colour = "black", fill = NA),
panel.background = element_blank(),
legend.position = "bottom",
plot.margin = unit(c(1.5, 1.5, 1.5 ,1.5 ),"cm")) +
facet_wrap(~data_type, ncol = 3, labeller = as_labeller(c("mumax" = "Growth rate",
"lambda" = "Lag phase",
"cv_area_batch" = "CV cell size",
"mean_AR_batch" = "Cell shape",
"mean_area_batch" = "Cell size",
"var_area_batch" = "Variance cell size")), scales = "free_y")
# print plot
response_plot
model_38_predictions %>% group_by(data_type) %>% select(-percent_batch1, -lineage) %>% filter(batches_38 == 3)
## # A tibble: 6 x 7
## # Groups: data_type [6]
## data_type data mixed_model batches_38 fit lwr upr
## <fct> <list> <list> <int> <dbl> <dbl> <dbl>
## 1 mumax <tibble [20 × 7… <MCMCglmm> 3 4.72 0.468 9.08
## 2 lambda <tibble [20 × 7… <MCMCglmm> 3 0.756 -0.614 2.15
## 3 cv_area_batch <tibble [20 × 7… <MCMCglmm> 3 32.9 5.75 56.6
## 4 mean_AR_batch <tibble [20 × 7… <MCMCglmm> 3 -16.5 -19.0 -14.1
## 5 mean_area_batch <tibble [20 × 7… <MCMCglmm> 3 -30.2 -49.8 -11.2
## 6 var_area_batch <tibble [20 × 7… <MCMCglmm> 3 -17.4 -39.5 4.79
model_38_predictions %>% group_by(data_type) %>% select(-percent_batch1, -lineage) %>% filter(batches_38 == 0)
## # A tibble: 6 x 7
## # Groups: data_type [6]
## data_type data mixed_model batches_38 fit lwr upr
## <fct> <list> <list> <int> <dbl> <dbl> <dbl>
## 1 mumax <tibble [20 × 7]> <MCMCglmm> 0 1.58 -2.30 5.68
## 2 lambda <tibble [20 × 7]> <MCMCglmm> 0 7.72 6.46 9.09
## 3 cv_area_batch <tibble [20 × 7]> <MCMCglmm> 0 43.5 19.6 69.2
## 4 mean_AR_batch <tibble [20 × 7]> <MCMCglmm> 0 -10.7 -12.9 -8.38
## 5 mean_area_batch <tibble [20 × 7]> <MCMCglmm> 0 -27.3 -47.5 -9.27
## 6 var_area_batch <tibble [20 × 7]> <MCMCglmm> 0 13.5 -6.14 36.0
model_20_predictions %>% group_by(data_type) %>% select(-percent_batch1, -lineage) %>% filter(batches_20 == 1)
## # A tibble: 6 x 7
## # Groups: data_type [6]
## data_type data mixed_model batches_20 fit lwr upr
## <fct> <list> <list> <dbl> <dbl> <dbl> <dbl>
## 1 mumax <tibble [8 × … <MCMCglmm> 1 0.442 -0.438 1.47
## 2 lambda <tibble [8 × … <MCMCglmm> 1 0.557 -0.607 1.64
## 3 cv_area_batch <tibble [8 × … <MCMCglmm> 1 -10.3 -13.6 -6.62
## 4 mean_AR_batch <tibble [8 × … <MCMCglmm> 1 -5.16 -10.6 0.0956
## 5 mean_area_batch <tibble [8 × … <MCMCglmm> 1 -11.1 -16.6 -5.01
## 6 var_area_batch <tibble [8 × … <MCMCglmm> 1 -36.9 -43.8 -30.5